% OLS and Newey West standard error
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [b tstat tstat_NW adjR2,res,DW,F]=ols(y,x,K,R)
% compute ols estimate and consistent standard errors
% of y = xb + u
% K is the autocorrelation order for NW
% K==0 equals white corrections
% R is the set of linear restrictions to test
% R=[0 0 1 0;
%    0 0 0 1];
% tests whether b_3 = b_4 = 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [b,tstat,tstat_NW,adjR2,res,DW,F] = OLS(y,x,K,R)

[n k]=size(x);
b=x\y;
fit=x*b;
res=y-fit;
W0=(res'*res)*inv(x'*x);
se=diag(sqrt(W0/(n-k)));
tstat=b./se;
%%%%%%%Newey West%%%%%%%
if K>=0
    xu=x.*repmat(res,1,k);
    kn=-K:K;
    V=zeros(k,k);
    for i=1:length(kn)
        Gamma=zeros(k,k);
        for j=max(1,kn(i)+1):min(n,n+kn(i))
            Gamma = Gamma+ xu(j,:)'*xu(j-kn(i),:);
        end
        Gamma = sum(Gamma,3);
        v=(1-abs(kn(i))/(K+1))/(n-kn(i))*Gamma;
        V=V+v;
    end
    W=inv(x'*x/n)*V*inv(x'*x/n);
    se_NW=diag(sqrt(W/n));
    tstat_NW=b./se_NW;
else
    tstat_NW=[];
    se_NW = [];
end
%%%%%%adjR2%%%%%%%%%%%%
R2=sum((fit-mean(y)).^2)/sum((y-mean(y)).^2);
adjR2=(1-k)/(n-k)+(n-1)/(n-k)*R2;
%%%%%Durbin-Watson%%%%%
rr=res(2:end)-res(1:end-1);
DW=rr'*rr/(res'*res);
%%%%%F-test on joint significance%%%%%
if exist('R','var')
    if ~exist('c','var')
        c=0;
    end
    F=n*(R*b-c)'*inv(R*W0*R')*(R*b-c);
    if F>chi2inv(.95,size(F,1)-1)
        disp('F-test Null is rejected')
    end
end

return